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Abstract 

O ! We study the discretization of the Escape Time problem: find the length of the shortest 

1} path joining an arbitrary point z of a domain Q, to the boundary dQ. Path length is 

measured locally via a Finsler metric, potentially asymmetric and strongly anisotropic. This 
Optimal Control problem can be reformulated as a static Hamilton Jacobi, or Anisotropic 
Eikonal, Partial Differential Equation, as well as a front propagation model. It has numerous 
applications, ranging from motion planning to image segmentation. 
^ We introduce a new algorithm, Fast Marching using Anisotropic Stencil Refinement (FM- 

ASR) , which addresses this problem on a two dimensional domain discretized on a cartesian 
grid. The local stencils used in our discretization are produced by arithmetic means, like 
in the FM-LBR [5], a method previously introduced by the author in the special case of 
Riemannian metrics. The complexity of the FM-ASR, in an average sense over all grid 
orientations, only depends (poly-) logarithmically on the anisotropy ratio of the metric, while 
i i most alternative approaches have a polynomial dependence. Numerical experiments show, 

in several occasions, that the accuracy /complexity compromise is improved by an order of 
magnitude or more. 
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Introduction 



qq The Escape Time D(z), from a point z of the domain $7, is the length of the shortest path joining 

this point to the boundary dd. Computing the escape time, and extracting an associated 
minimal path, is a task of obvious interest in motion planning control problems [12]. Yet 
• • this versatile problem has numerous other applications |15j . including image classification |17j . 

. £h seismic imaging [5] or the modeling of bio physical phenomena |16| . We are motivated by medical 

image segmentation problems, which often involve a strongly anisotropic [TJ, and potentially 
^ asymmetric |X3|, I14j . local measure of path length. 

From a theoretical point of view, the Escape Time problem can be reformulated as a static 
Hamilton Jacobi, or Anisotropic Eikonal, Partial Differential Equation (PDE) [15]. Its numerical 
discretization has attracted a lot of research effort, and includes the Fast Marching algorithm |10| . 
the Fast Sweeping method [18], and their numerous variants [5j [T2], [2] . As the "Fast" adjective 
indicates, performance is a crucial concern: in image processing applications, the discretization 
domain may contain millions of points (as many as image pixels), and CPU time should remain 
compatible with user interaction. Last but not least, as mentioned above, state of the art image 
processing applications involve strongly non- uniform, anisotropic and/or asymmetric measures 
of path length, which challenges available algorithms pQ and limits the parallelization potential 
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This paper is devoted the introduction and study of a new algorithm, Fast Marching using 
Anisotropic Stencil Refinement (FM-ASR), a numerical solver for the two dimensional Escape 
Time problem discretized on a cartesian grid. Path length is measured locally through a given 
arbitrary Finsler metric F: a continuous map associating to each point z G fl an asymmetric 
norm F z . The FM-ASR regards the discretization grid as a subset of the Lattice 7L 2 , and 
uses arithmetic tools to produce the local stencils involved in the discretization of the associated 
Partial Differential Equation (PDE), which results in a huge complexity reduction in comparison 
with more classical approaches. Note that the FM-LBR |5j, previously introduced by the author, 
shares this approach but is limited to metrics of Riemannian type (elliptic anisotropy). Non- 
Riemannian metrics arise for instance in applications which take advantage of their potential 
asymmetry [141 [13] , or as the result of the homogenization of smaller scale Riemannian metrics 
[7]. The anisotropy ratios of an asymmetric norm F : H 2 — > H + , and of a Finsler metric 
F : H x R 2 -> R + , are defined by 

F(u) 

k(F) := max — -, k{F) := sup k(F z ). (1) 
H=|«|=i b [v) zeU 

The average complexity of the FM-ASR only depends (poly-) logarithmically on the anisotropy 
ratio of the given metric F, and is quasi-linear in the number N of discretization points. In 
contrast, alternative approaches show a polynomial dependence either on k(F) |12[ [H], or on 
N [2], a difference clearly apparent in the numerical experiments presented in §3. In average 
over all grid orientations, and denoting ln a x := (lnx) a , the complexity of the FM-ASR is only 
0{N\i^ k{F)+ N\nN). 



1 Description of the problem, algorithm, and main results 

The Escape Time problem is posed on a two dimensional bounded domain ft C K 2 , equipped 
with a Finsler metric F. This metric is a continuous map F : Q x H 2 — > H + , (z, u) i— > J- Z (u), 
such that for each fixed z G O, the restriction u \— > F z (u) is an asymmetric norm (i.e. a proper 
homogeneous convex function^]) . The length of a path 7 G C 1 ([0, 1], Q) is measured through the 
metric T: 

length( 7 ) := / F l[t) {i {t)) dt. 
Jo 

Notable special cases include Isotropic metrics: F z (u) = n(z)\\u\\, where the parameter n(z) > 
corresponds to the local index in geometrical optics. Riemannian metrics have the form: 
F z (u) := \J (u, J\4(z)u) , where A4(z) is a symmetric positive definite matrix. Symmetric Finsler 
metrics are subject to the condition F z (—u) = F z (u), for all z G ft, u G K 2 . See Figure [l] Here 
and below we denote by || • || and (•, •) the canonical euclidean norm and scalar product on H 2 . 

The length of a path 7 G C 1 ([0, and of the reversed path 7 : t 1— > 7(1 — t) may be 

different in the case of a general asymmetric Finsler metric. This apparent oddity is entirely 
relevant in the study of motion planning under the influence of wind |12j . It is also essential 
in minimal path based image segmentation methods |144 113] . where the right and left of the 
path should have different prescribed characteristics, since they respectively correspond to the 
foreground and background of the segmented object. We introduce an asymmetric distance 
D(-,-) 011O 

D(x,y) := inf{length( 7 ); 7 G C\[0, 7 (0) = x, 7 (1) = y}. 

1 Finsler metrics are often assumed to be smooth, and the local asymmetric norms to be strictly convex. These 
assumptions, tailored for the study of minimal paths, are not required in our analysis of the Escape Time problem. 



2 






Figure 1: A Finsler metric J 7 , on a domain C H 2 , is the data of a continuously varying 
asymmetric norm F z , at each point z £ Q. The convex sets {u; J- Z {v) < 1}, at several points 
z G $7, are used to visualize the metric T . Finsler metrics can be of Riemannian type (left), 
symmetric (center), or asymmetric (right). The discretization of the Escape Time problem 
involves the construction of local stencils, three those produced by the FM-ASR are illustrated. 



The solution of the Escape Time optimal control problem is the distance D(-) to the boundary: 
for all x G O 

D(x) := min{D(x, y); y G dtl}. (2) 

The function D is also characterized as the unique viscosity solution [3] of the static Hamilton 
Jacobi, or Anisotropic Eikonal, PDE (see e.g. [7] for a discussion of this reformulation) 



J"*(-V D(z)) = 1 for all ze!l, 
B(z) = for all z G <9R 



(3) 



In the above equation, we denoted by F* the dual asymmetric norm of an asymmetric norm F 
on 1R 2 , which is defined for all u G H 2 by 



F*(u) := max 



(u,v) 



v^o F(v) ' 



(4) 



Consider the front defined by £t := -D -1 ({£}), t > 0, thus So = <9fl The normal to this front, 
at a point z G £t where D is differentiable, is positively collinear to VD(z). The speed of the 
front along in this normal direction is inversely proportional to the gradient euclidean norm, 
l/||VD(z)||, and is thus determined by the identity J'*(-VD(z)) = 1. Note that the front may 
only go forward, and that the front speed cannot depend on global or high order properties of 
the front, such as its curvature. See [15] for the applications, and limits, of this elementary front 
propagation model. 



Since D(-, •) is a path length (asymmetric) distance, one has for any point x and neighborhood 
V, x G V C O, the identity 

D(x) = min D(x, y) + D(y). (5) 

y£dV 

Indeed, any path 7 joining x to d£l must cross dV at least once, at some point y. The discretiza- 
tion of the Escape time problem is based on an approximation of the right hand side of ([5]), the 
so-called Hopf-Lax update operator introduced in [20] , see also [HJ EJ [5] , and a reinterpretation 
of this equation as a fixed point problem. 

For that purpose we introduce discrete sets f2* and <9f2*, devoted to the sampling of the 
continuous domain Q and of its boundary d£l respectively. In the FM-ASR, and d£l* need 
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Figure 2: Stencil used in the classical Fast Marching algorithm (first), the AGSI [2j (second), 
the FM-8 (third), and the MAOUM [12J at a grid point z € f2* when the local anisotropy ratio 
n{J- z ) is 1.5 or 6 (fourth and fifth, respectively). Algorithms compared in ^4] to the FM-ASR. 



to be subsets of the grid 2Z 2 , or of another orthogonal grid obtained by rescaling, rotating and 
offsetting 2Z 2 . A small neighborhood V#(z) of each z £ !!„ the stencil, is constructed under 
the form of a triangulation, of vertices in U d£l* . See Figure [T] for some stencils used in the 
FM-ASR, and Figure [2] more classical example^] For any discretization point x G f2*, and any 
discrete map d : f2* U — > 1R+, we define the Hopf-Lax update 

A(d,a;) := min F x (y - x) + l v d(y). (6) 

We denoted by V the stencil V*(x), and by Iy the piecewise linear interpolation operator on 
this triangulation. Note that A(d, x) does not depend on the value of d(x), but only on d(y) for 
points y of the discrete domain 0* U <9f2* which lie on the boundary of the stencil V = V*(x). In 
the following we set 1R + := IR+ U {+00} = [0, +00], adopt the convention x 00 = 0, and allow 
discrete maps to take the value +00. 

Numerical methods for the Escape Time problem construct a discrete approximation d : 
SI* U — > 1R+ of the continuous solution D of ^ , characterized by the following discrete fixed 
point problem: 

d(z) = A(d, z) for all z£l] t , , . 

d(z) = for all dO*. ^ ' 

Note that the distance on a weighted graph obeys a system of equations of similar nature, 
except that the neighborhoods V*(z), of each vertex z, are given by the graph structure and 
not dependent on the numerical method. Two well known algorithms can be used to solve this 
system and evaluate graph distances: the fast, single pass, Dijkstra algorithm, and the slower 
but more flexible (in that negative edge weights are allowed) Bellman-Ford algorithm. The 
algorithms used to solve the system (J7|), associated to the Escape Time problem, are inspired 
by these two methods, and the lack of negative edge weights in the graph setting is translated 
into a geometrical property of the stencils, named the Causality Property, see below. 

Bellman-Ford inspired algorithms solve the system Q via Gauss-Siedel iteration: the replace- 
ment rule d(2%) ::= A(d, z^), k > 0, is applied repeatedly to a mutable map d : fi* U dQ* — > 1R+, 
until a prescribed convergence criterion is met. The map d is initialized to +00 on fi*, and 
on 8ft*. The choice of the sequence of points z^ 6 k > 0, depends on the method. This 
sequence enumerates the lines and columns of in the fast sweeping methods [18] , and is ob- 
tained via a priority queue in the Adaptive Gauss Siedel Iteration (AGSI) [2j. The stencils are 
usually extremely simple, see Figure [2j left and center left. The complexity of these methods is 
linear in N := #(^2*) in the special case of an Isotropic metric, 0(\(J-)N) for the Fast Sweeping 

2 The stencil construction of the AGSI and of the MAOUM requires a mesh of the underlying discrete domain 
fi*, here a subset of hTZ 2 for some h > 0. We triangulated this grid with rescaled translates of the triangle of 
vertices (0, 0), (1, 0), (0, 1), and of its symmetric with respect to the origin. 



[22], but is polynomial in general, 0(/i(J")iV 3 / 2 ) for the AGSI [2]. The constants A (J 7 ) and fj,(T) 
depend on global geometrical features of the metric. The AGSI is popular, simple and quite 
efficient; it appears for reference in our numerical experiments. 

We next introduce some geometrical concepts, and the Causality Property which is at the 
foundation of Dijkstra inspired solvers of the Escape Time problem: the Fast-Marching algorithm 
[TO] , and its variants [HI [T2l [5] . When satisfied, this property allows to "decouple" and solve the 
discrete system (j7| in a non-iterative, single pass fashion, resulting in a complexity independent 
of global features of the metric, and quasi-linear in the number N of unknowns. 

Definition 1.1. Let F be an asymmetric norm on M 2 . We say that two vectors u, v G M 2 \ {0} 
form an F -acute angle if 

F(u + 5v) > F(u) and F(v + 5u) > F(v) for all 5>0. (8) 

We say that a finite conforming triangulation T is F -acute if 

(i) The union of the triangles T £ T is a neighborhood of the origin. 

(ii) The vertices of each T G T lie on 7L 2 , one of them is the origin 0, and T has area 1/2. 
(Hi) The non-zero vertices of each triangle T G T form an F -acute angle. 

In other words, two vectors form a F-acute angle if adding a positive multiple of one to the 
other increases its F norm. Condition (i) heuristically ensures that the information is propagated 
in all directions in 0. Condition (ii) ensures that this information stays on the grid £7*. In 
addition, this condition implies that the union of the triangles T G T does not contain any point 
of 7L 2 in its interior except the origin, which heuristically ensures that information does not 
"jump over" a subset of fi*. 

The core of this paper is devoted to the construction and study of a -F-acute mesh T(F), 
defined for each asymmetric norm F, see Figure [TJ [3j and used to assemble the stencils of the 
FM-ASR. This mesh is produced by the following algorithm. 

Construction of the mesh T(F), associated to a given asymmetric norm F. 
This mesh is star shaped with respect to the origin, see Figure [TJ The sequence L of its 
consecutive boundary vertices is generated as follows, using only two lists L and M. 
Set L :: = [(1, 0)], M ::= [(1, 0), (0, -1), (-1, 0), (0, 1)]. 
While M is non-empty do 

Denote by u, v the last element respectively of L and M. 
If u, v form a i^-acute angle 

then remove v from M and append it to L 
else append u + v to M. 



Assume that the discrete domain f2* is a portion of the grid Z := zq + hRTL 2 , obtained by 
rotating the grid 2Z 2 by a rotation R, rescaling it by a factor h > 0, and offsetting it by zq G H 2 . 
Then the stencil V*(z), assembled in the Preprocessing of the FM-ASR, algorithm page 

[6j and involved in ([7]), is defined by rotating, rescaling and offsetting the mesh T(F Z o R): with 
obvious notations 

V*(z) :=z + hRT{F z oR). (9) 

These stencils have a fine angular resolution in the direction of anisotropy, and a coarser one 
in other directions, see Figure [TJ This distinctive property justifies the name of our algorithm: 
Fast Marching using Anisotropic Stencil Refinement (FM-ASR). 
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Figure 3: The mesh T{F) constructed for norms F of anisotropic euclidean type (top, i.e. given 
by a matrix M G S^), of anisotropy ratio k(F) ranging from 1 to 32 and orientation tt/3 (top 
left), or anisotropy ratio k(F) = 10 and orientation ranging from 7r/4 to tt/2 (right). Likewise 
for asymmetric norms of the form F(u) := \\u\\ — (oj,u), u G H 2 , (bottom) of anisotropy ratio 
ranging from 4 to 400 and orientation tt/3 (bottom left), or of anisotropy ratio 100 and varying 
orientations (bottom right). 



FM-ASR: Preprocessing. 
Input: A bounded domain £7 C H 2 , equipped with a Finsler Metric J £ C°(S) x H 2 ,1R + ). 

A grid Z, obtained by rotating, rescaling and offsetting (if needed) the grid 2Z 2 . 
Set a* :=fin2. 

Assemble the stencils V*(z), z G $7*, as in 

Assemble the "reversed stencils", defined by V*(y) := {x G r2* \ {y}; y is a vertex of V*(x)}. 
Setdn*:={yeZ\n*;V*(y)^®}. 



FM-ASR: Execution (common to other variants of the Fast Marching algorithm |10j). 
Variables: a boolean table b : SI* U — > {trial, accepted}, and a map d : SI* U 517* — > H + . 
Initialize d to +oo on SI*, and to on 3S1*. Initialize b identically to trial. 
While b is not identically accepted do 

Denote by y G S7* U <9S1* a minimizer of d among those points such that b(y) = trial. 
Set b(y) ::= accepted. 

For all x G V*(y) such that b(y) = trial do 
Set d(x) ::= min{d(a;), A(d, x; b,y)}. 
Output: the distance map d. 
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The vertices v of the mesh T(F), associated to an asymmetric norm F, are bounded in terms 
of the anisotropy ratio: ||t> || < 2k(F) (see Proposition |2.9| below) . Hence the discrete boundary 
<9f2*, produced by the FM-ASR initialization, may contain grid points at distance 2/ik(J 7 ) from 
the domain Q. This is not an issue in the case of the null boundary condition ([3]), Q, chosen in 
our presentation, or of a point source problem (the most common case in applications, see 
However if the boundary condition is non-trivial, then its extension from the boundary d£l to 
the wider discrete set dQ* is required by the FM-ASR. 

The next lemma gives a simple characterization of -F-acuteness when the asymmetric norm 
F is differentiable or of anisotropic euclidean type (i.e. defined by a symmetric positive definite 



matrix). The characterization (10), for smooth norms, was introduced in [21J in the same 
context. We denote by the collection of 2 x 2 symmetric positive definite matrices. 

Lemma 1.2. Let F be an asymmetric norm on M 2 , and let u,v G R 2 \ {0}. 

1. If F is differentiable at u,v, then these vectors form an F -acute angle if and only if 

(u,VF(v))>0 and (v,VF(u))>0. (10) 

2. If these exists M G S^" such that F(w) = \f (w, Mw), for all w G M 2 , then u,v form an 
F -acute angle if and only if 

(u,Mv)>0. (11) 

Proof. We first establish Point 1. We have the Taylor development F(u + 5v) = F(u) + 
5(v,V F(u)) + o(5) as 5 — >• 0, and likewise exchanging the roles of u and v. Thus ^ clearly 



implies (10). Conversely the function F, being convex, is above its tangent maps, hence 



F(u + 5v) > F(u) + 6(v, VF(u)) for all S G IR, and likewise exchanging the roles of u and 



v. Thus (10) implies (|8j), which concludes the proof of Point 1. 

Point 2 immediately follows from the following expansion: for any u, v G TR 2 , 5 G H, one has 

F(u + 5v) 2 = F(u) 2 + 26{u, Mv) + 5 2 F{v) 2 . □ 

If F is the canonical euclidean norm, then F-acuteness coincides with the standard notion 



of acuteness (apply (11) to M := Id). The following proposition, or a close variant (HUE], is 
at the foundation of all Dijkstra inspired methods. The positivity of the differences d w — d u , 
d w — d v , is a substitute for the positivity of the edge weights in the classical Dijkstra algorithm. 

Proposition 1.3 (Causality Property). Let F be an asymmetric norm on M 2 , let u,v G M 2 be 

linearly independent, and let d u ,d v G M. Assume that u and v form an F -acute angle. Define 

d w := min td u + (1 - t)d v + F(tu + (1 - t)v), (12) 
te[o,i] 

and assume that this minimum is not attained for t G {0, 1}. Then d u < d w and d v < d w . 

Proof. See appendix. □ 

In order to describe the Execution of the FM-ASR, see algorithm page [6j we introduce 
a variant of the Hopf-Lax update ([6]), which uses two additional variables: a boolean map 
b : (7* — > {trial, accepted}, and a grid point t/ 6 Sl t . 

A(d, x; b, y) := min T{x — x) + ly d(x'), (13) 
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where V := V*(x), and T denotes union of the vertex y, and the (at most two) segments [y, z] of 
dV containing y and another vertex z of V such that b(z) = accepted. This second part of the 
FM- ASR is common to the original Fast-Marching algorithm |10j and its variants |12|. [5] . The 
fact that it solves the discrete fixed point problem ([7]) follows from the Causality Property, see 
[TUl [5] for a proof. 

For each step size h > 0, consider the discrete domain $7^ := n (h7L 2 ), and the associated 
solution dh of the system Q produced by the FM-ASR. A proof of uniform convergence of the 
discrete maps (dh)h>o towards the solution D of the continuous Escape Time problem, 



lim ( max I D(z) — dh(z)\ } 



h 

was obtained in [5] for the FM-LBR, a closely related algorithm, in the special case where O is 
the unit square equipped with periodic boundary conditions, and dVt is reduced to a point. The 
adaptation of this proof to the FM-ASR is straightforward^ and is not reproduced here. 

We next turn to the analysis of the complexity of the FM-ASR, beginning with the Pre- 
processing step. We do not count in the complexity the construction of the discrete domain 
O*, as the intersection of a grid Z with the continuous domain f2, since this is either trivial 
or dependent on the chosen representation of f2. Consider an asymmetric norm F such that 
answering the predicate "it, v form a .F-acute angle" has cost 0(1), for any u, v G 1R 2 . That is 
the case is F is differentiable, and if evaluating the gradient VF(u) has cost 0(1), using Lemma 



1.2 Then constructing the -F-acute mesh T(F) has cost 0(#(T(F))), where #(T(F)) denotes 
the number of triangles in the triangulation T(F), which is also the number of its boundary 
vertices. As a result, assembling the stencils of the FM-ASR has cost O(N'), where 

Assembling the reversed stencils V*(z), z £ £1*, is done by reversing a directed graph having N' 
edges, and thus also has cost O(N'). Note that N' is also the sum of the cardinalities of the 
reversed stencils. Storing these stencils leads to a O(N') memory footprint for the FM-ASR, 
which is not required by e.g. the AGSI |2j, see Remark 2.5 in [5 J for a discussion of this point. 
The total complexity of the FM-ASR Preprocessing is thus O(N'). 

We next turn to the FM-ASR execution, and for that purpose we denote by N := $:(Sl*U<9f2*) 
the cardinality of the discrete domain, which is also the number of unknowns in ([7]). The ex- 
ecution requires to maintain a list of the points in fi* U <90* such that b(z) = trial, sorted by 
increasing values of d. We assume in this complexity analysis that the data structure used for 
this purpose is a Fibonacci Heap, in such way that the "Remove -Key" and "Decrease_Key" oper- 
ations on this list have respective amortized complexity 0(ln N) and 0(1). The "Remove -Key" 
routine is called N times, once a each command b(y) ::= accepted, and the "DecreaseJKey" 
routine at most N' time^J once at each command d(x) ::= min{d(x), A(d, x; b,y)}. Evaluating 
the modified Hopf-Lax update operator A(d, x; b, y) requires to solve at most two convex min- 



imization problems of the form (12): one for each boundary edge [y,z] of V*(x) containing y 



and a vertex z such that b(z) = accepted (13). The complexity of their resolution is regarded as 



3 The proof can in fact be simplified in the case of the FM-ASR, since the stencil V*(z) of a grid point z £ fi* 
contains the four immediate grid neighbors of z. This makes Lemmas 2.6 (Consistency) and 2.7 trivial in 0. 

4 Fibonacci Heaps are a data structure specifically tailored for Dijkstra-Like algorithms on densely connected 
graphs: N' 2> N. In the numerical experiments presented on Q one always have N' < 20iV for the FM-ASR, and 
using a classical binary heap proved to be more efficient. We used Boost's implementation of Fibonacci heaps, 
and the Standard Template Library for binary heaps. 
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elementary; in many interesting cases, they have an explicit solution involving 0(1) elementary 



operations (+,— , X,/ and among reals, see Proposition 4.1. Like other Dijkstra inspired 
algorithms, the total cost of the FM-ASR execution is thus 0(N' + NlnN). 

The rest of this paper is devoted to the estimate of total the stencil cardinality N' of the 
FM-ASR, or more precisely of its average value jV' over all rotations and translations of the 
discretization grid. We show that iV' < N(l + ln a n(F)) with a = 2 for symmetric Finsler 
metrics, and a = 3 otherwise, see Figure [l| The average case complexity of the FM-ASR is 
thus 0{Nln a n{F) + NlnN). This is significantly below Bellman-Ford inspired algorithms, 
such as the AGSI 0{\{F)N 3 / 2 ) [2], thanks to the quasi-linear complexity in N. Alternative 
Dijkstra inspired solvers include the Ordered Upwind Method (OUM) [8] (which uses dynamic 
stencils, constructed on the fly during the execution), and the Monotone Acceptance OUM 
(MAOUM) |12j . They use larger stencils, of cardinality between k(F) and k(F) 2 , which results 
in a complexity linear if not polynomial in the anisotropy ratio: 0(k(F)^N In N) [U [12], for 
som^] {3 G [1,2]. Fast Marching using Lattice Basis Reduction (FM-LBR), introduced in [5] 
by the author, has like the FM-ASR a complexity 0{N In k{F) + NlnN) logarithmic in the 
metric anisotropy and quasi linear in the number of unknowns. Yet the application range of 
the FM-LBR is different: it extends to dimension 3, and 4 [6], but only applies to metrics of 
Riemannian type. A numerical comparison of the FM-ASR with the AGSI, MAOUM, FM-8 
(a fast but not always convergent alternative) and FM-LBR (when applicable) , is presented in Q 

The following theorem is the main result of this paper. It establishes that the cardinality of 
the -F-acute mesh #F(F), used to define the stencils of the FM-ASR, grows only logarithmically 
with the anisotropy of the given asymmetric norm F, in an average sense over all orientations. 
For each 9 G H we define the rotated asymmetric norm F e by 

F e {u) := F(Rju), 

where u G 1R 2 and Rq denotes the rotation matrix of angle 9, see Figure |4j 

Theorem 1.4. There exists a constant C, such that for any asymmetric norm F on M 2 , one 
has: 

i-2-k 

/ #(T(F 9 ))d9 <C(1 + In 3 k(F)) (14) 
J o 

A slightly sharper estimate holds if F is symmetric: 

i-2-k 

/ #(T(F e ))d9 <C(1 + In 2 k(F)). (15) 
J o 

This result immediately translates into an estimate of the average complexity of the FM- 
ASR. Indeed consider a fixed discretization step h > 0. For all 9 G [0, 2tt[ and all u G [0, 1[ 2 , let 
Zq^ u be the image of the grid 7L 2 by the affine map : z ^ hRj (z+u). Let Qg iU := Q,C\Zq u be 
the discrete domain, and let dVLg^ u be its discrete boundary, see the FM-ASR Preprocessing page 
pi The number of discretization points #($7e jU U dQg jU ) is close to iV := |f2|/i~ 2 (independently 



5 Strictly speaking, /3 = 2. Yet the asymptotic complexity of the OUM, as N — ¥ oo, drops to 0(n(J r )N \n N). 
If the MAOUM is executed on a periodic mesh, then the stencil of z only depends on a single parameter: 
the anisotropy ratio k(T z ). It has cardinality 0(K(J r z )), but costs 0(k(J-" z ) 2 ) to construct. In our numerical 
experiments ij4] these stencils are precomputed, stored in a look-up table, and the complexity of the MAOUM 
drops to 0{n{T)N In N). 
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Figure 4: The unit ball of F is the unit ball of F rotated by the angle 9 (left). Here the norm 
F, of anisotropic euclidean type, is given by the diagonal matrix of entries (k, 1/k), with k = 4 
(left), k = 100 (center, linear plot) and k = e 8 (right, log plot). In this example, the cardinality 
of T{F e ) is highly dependent on the angle 9, and seems to spike when (cos 9, sin 9) is close to 
be proportional to a vector with small integer coordinates. 



of 9 and u), provided the domain f2 is smooth and h is sufficiently small. The average value N' 
of the sum of the stencil cardinalities, as 9 £ [0, 2ir] and u G [0, 1[ 2 , is given by 

r2n p pin r r 

/ / J2 #(F 9 z )dud9 = h- 2 / #(T(F 9 z ))dzd9 < Ch~ 2 (1 + ln a n(F z ))dz, 
Jo J[o,i[ 2 Jo Jzen Jzen 



where the constant C and the integer a G {2,3} are as in Theorem 1.4 The first equality is 
obtained by a change of variables, observing that the affine transformation ipg tU has Jacobian 
determinant h 2 , and that the sets (^e,n) u e[o,i[ 2 partition £2, for any fixed 9 G [0, 2tt]. Hence 
N' < CN{1 + ln a k(F)) as announced. ' 

The next proposition discusses worst case performance. In this case, the impact of anisotropy 
on mesh cardinality can be quasi linear, instead of logarithmic, and our anisotropic construction 
of F-acute meshes T(F) has little advantage over an isotropic one T k (f)-, depending only on the 
anisotropy ratio. This proposition implies the worst case upper bound N' < Nk(J-)(1 + \u.k(J-)) 
on the total stencil cardinality, N' < Nk(F) for Riemannian metrics. The worst case complexity 
of the FM-ASR is thus 0(Nk(F) In k(F) + N\nN), or 0(Nk(F) + N\nN) for Riemannian 
metrics, which is comparable to the OUM and MAOUM. 

Proposition 1.5. There exists a constant C, such that for any asymmetric norm F on M 2 one 
has: 

#(T(F))<CK(F)(l + ln K (F)). (16) 
A slightly sharper estimate holds if F is symmetric: 

#(T(F)) < Ck(F). (17) 

For any k > 1, there exists a "non- adaptive" mesh T K which is F-acute for any asymmetric 
norm such that n(F) < k, and has cardinality #(7^) < Ck{1 +lnre). There also exists an 
anisotropic euclidean norm F K such that k(F) < n and #(T(i ? K )) > n/C. 

We discuss in ^2] the construction of the i^-acute mesh T(F), for any asymmetric norm F. 
Section §3 is devoted to the proof of our main result Theorem |1.4[ achieved in Corollaries 3.10 
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and 3.13 The proof of the worst case analysis, Proposition 1.5 is achieved in Corollaries 2.10 
and 3.10 We present some numerical results in Q 

Remark 1.6 (Performance comes at the price of specialization). The FM-ASR, introduced in 
this paper, is an efficient method to solve strongly anisotropic and/or asymmetric Escape Time 
problems when the discrete domain f2* is a subset ofTZ? , or of another orthogonal grid. Extending 
this algorithm to a broader class of discrete domains requires to generalize the construction of 
the stencils V(z), z G fi*. In particular one must find an analog of the rule "ifu,v G T2? do not 
form an T z -acute angle, then consider their sum u + v G 7Z?" which appears implicitly in the 
construction of the mesh T(J- Z ), algorithm page^ and thus ofV(z) This is non-trivial. 

• // the discrete domain f2*, two dimensional, is not a grid subset. The points u,v 
do not belong to a lattice, but are differences u = x — z, v = y — z, between the point z where 
the stencil is constructed, and close-by discrete points x,y G f2*. The new inserted stencil 
vertex z' cannot be obtained as the sum z + u + v = x + y — z, which may not belong to . 
Instead, z' should be chosen as the point closest to z in the open cone z + M+u + IR+v, 
where IR* + denotes positive reals. However, it is not clear wether data structures exist, for 
the discrete domain 0,*, which allow to perform this closest point search without strongly 
increasing the complexity of the FM-ASR. 

• // the domain is three dimensional, and 17* a subset of H? . There are now three 
points u,v,w G Tl? , vertices of a facet of the stencil boundary dV. The extension of 
the Causality Property, Proposition \1.3\ in the case of an arbitrary asymmetric norm 
F : IP? — > M+, requires not only the pairs of vertices (u, v), (u, w), (v, w) to form F -acute 
angles, but also all the pairs of a vertex and a point of the opposite edge: (u, tv + (l — t)w), 
t G [0,1], and likewise exchanging the roles of u,v,w. This can be be difficult to check 
numerically. In the case of a Riemannian metric, checking F-acuteness for pairs of vertices 
is sufficient fS\, but there remains an ambiguity: should the new inserted vertex beu + v 
or u + w, if none of the corresponding angles is F -acute? Our attempts to generalize the 
FM-ASR to this setting were unconvincing, both experimentally and theoretically, hence 
we recommend the FM-LBR [5] for such 3d, Riemannian, Escape Time problems. 

2 Construction of the stencils 

We discuss in this section the construction of the F- reduced mesh T(F) , defined for each asym- 
metric norm F, and used to define the stencils of the FM-ASR ([9]). The construction presented 
in the introduction is reformulated as an in-order transversal of four binary trees. We establish 



a worst case upper bound on ff(T(F)) in Corollary 2.10 and we introduce a number of tools 
that will be used in ^3] to estimate the average cardinality of T(F e ), 9 G [0, 2-ir]. 

2.1 Mesh generation by recursive refinement 

All the triangles considered in the rest of this paper share some properties of geometric nature (or 
arithmetic nature, depending on the point of view), which are introduced in the next definition. 

Definition 2.1. An elementary triangle T, is a triangle satisfying the following properties: 

• One of the vertices of T is the origin, and the the other two belong to 7L 2 . 

• Denoting by u, v the non-zero vertices of T, one has 

\det(u,v)\ = 1 and s(T) := (u,v) > 0. (18) 
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Figure 5: Refinement of a triangle (left). Mesh 7o (center left). First levels of the binary tree 
(center right) defined by the recursive refinements of T\. Mesh defined by the ASC: a s(T) > 5" 
(right). 

The second point of this definition can be rephrased as a geometrical statement: T has area 
1/2, and has an acute angle at the origin. Let us recall that for any two vectors u,v G R 2 one 
has the identity 

(u,v) 2 + det(u,v) 2 = |H| 2 |H| 2 . (19) 

If u, v are non-zero, and if it and —v are not positively collinear, we denote by < (it, v) G (— ir, it) 
their oriented angle: 

cos(<(«,«)) = ,, , „ ,, and sin(<(it, v)) = „ ,„ ./ . (20 

1 1 1£ 1 1 1 1 f 1 1 

The scalar product s(T) associated to an elementary triangle T reflects its thinness, indeed if 



u,v are its non-zero vertices then combining (18), (19) and (20) we obtain 



sin|<(7x,u)| = „ „ = (a(T) 2 + l)-i (21) 

We introduce in the next definition the refinement of an elementary triangle T, which is 
illustrated on Figure [5] (left). Note that T is strictly covered by the union of its children. 

Definition 2.2. The refinement of an elementary triangle T of non-zero vertices u,v consists of 
the two elementary triangles T' and T" of non-zero vertices (it, u+v), and (u+v, v), respectively, 
which are referred to as its children. 

The scalar product s(-) grows with refinement: 

s(T') = ( u ,u + v) = S{T) + ||it|| 2 > S(T) + 1. (22) 



This property, combined with (21) reflects the fact that the recursive children of an elementary 
triangle become thinner an thinner, as can be observed on Figure [5] (center right). 

We denote by 7o the mesh, illustrated on Figure [5j containing the four elementary triangles of 
non-zero vertices (±1,0) and (0, ±1). The next lemma establishes that any elementary triangle 
can be generated by recursive bisections from an element of 7o- 

Lemma 2.3. • LetT be an elementary triangle, of non-zero vertices u and v. The following 
are equivalent : (i) T G 7o, (H) \\u\\ = ||i>||, (Hi) (u,v) < min{||u|| 2 , ||i>|| 2 }. 

• The collection of elementary triangles, equipped with parent- children relationship, is a for- 
est of four infinite binary trees, which roots are the elements of To . 
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Proof. We begin with the first point. We clearly have (i) (££), by inspection of the four 
elements of 7o, and (ii) =>■ (m), by observing that it and v are not collinear and thus that 
(it, w) < ||u||||w||. We now assume (Hi) and establish (i). We have 

(u,u) 2 < min{||u|| 2 , |M| 2 } 2 < ||u|| 2 |H| 2 = (ii,v) 2 + det(u,u) 2 = (u,u) 2 + 1. 

Comparing the left and right and side, and observing that the members of these inequalities 
are all integers, we obtain that the large inequality above is an equality. Hence ||it[| = \\v\\ and 
(||n|| 2 ) 2 = (u,v) 2 + 1. Therefore (u,v) 2 and (||u|| 2 ) 2 are consecutive integers which are both 
perfect squares. Only the integers and 1 satisfy this property, hence 1 = ||u|| = \\v\\, which 
implies that these vectors are of the form (±1, 0) or (0, ±1). Since | det(it, v)\ = 1, these vectors 
are not collinear, and we obtain that T G 7o, which concludes the proof of the first point of this 
lemma. 



We next turn to the proof of the second point of this lemma. It follows from (22) that a 
triangle T 6 To cannot have a parent R, since it would satisfy s(R) < s(T) = 0. More generally, 
and for the same reason, an elementary triangle T has at most s(T) ancestors. 

In order to conclude the proof, we need to show that any elementary triangle T, which is 
not in 7o, has exactly one parent R. Let u, v be the non-zero vertices of T, ordered in such 
way that ||u|| < ||u||. The non-zero vertices of a parent R are either (it — v,v) or (u,v — u), 
but the first case can be excluded since (it — v,v) = (u,v) — \\v\\ 2 < \\u\\\\v\\ — |M| 2 < 0. 
Conversely, the triangle R which vertices are the origin, u and v — u is an elementary triangle 
since det(u, v — u) = det(n, v) = ±1, and (u, v — u) = (u, v) — \\u\\ 2 > 0. □ 

If a mesh T contains only elementary triangles, then it automatically satisfies assumption 
(ii) of Definition and so does any mesh T' obtained by refining, possibly recursively, some 
elements of T ■ If the mesh T satisfies assumption (i) of Definition |1.1| namely that the union of 
its elements is a neighborhood of the origin, then so does T' . The mesh constructions proposed 
in this paper consist in recursively refining the elements of a mesh until all of them satisfy a 
prescribed stopping criterion, see Figures [6] and [7J 

Definition 2.4. An admissible stopping criterion (ASC) is a predicate p which associates to 
each elementary triangle T a boolean value p(T), and which satisfies the following properties: 



(Heredity) LetT',T" be the children of an elementary triangle T. Ifp(T) holds, thenp(T') 
and p(T") also hold. 

(Finiteness) There exists a constant s p > such thatp(T) holds for any elementary triangle 
T satisfying s(T) > s p . 



The conjunction p A p' and the disjunction p\l p' of two ASCs p,p' are clearly also ASCs. 
We write p =>• p' if p(T) =>■ p'(T) for any elementary triangle T. 

Definition 2.5. Let p be an ASC. We denote by T(p) the collection of triangles obtained by 
recursively refining (i.e. replacing with their children) the elements of To, until each satisfies the 
predicate p. We denote by £(p) the collection of all elementary triangles which do not satisfy p. 



Definition 2.4 of an ASC p is tailored so that the set £ (p) can be identified to four finite 



binary trees within the forest of elementary triangles evoked in Lemma 2.3. The triangulation 
T(p) consists of the leaves of these trees. This is the main ingredient in the proof of following 
proposition. 
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Proposition 2.6. 1. The recursive procedure described in Definition \2.5\ ends after a finite 
number of steps, and yields a mesh T(p) which satisfies assumptions (i) and (ii) of Defi- 
nition \n\. Furthermore 

#(T(p)) = 4 + #(£(?)). (23) 

2. If two ASCs p,p' are such that p =>• p' , then #(T(p)) > #(T(p'))- 

3. For any two ASCsp,p', one has #(T(pAp')) < #(T(p)) + #(T(p'))- 

Proof. We denote by (Ti)i<i<4 the four elements of the mesh 7o, and by ("Pi)i<i<4 the four 
infinite binary trees evoked in the second point of Lemma 2.3, see also Figure [5} The root of 



Vi is the triangle T«, for any 1 < i < 4, and the children of any T E Vi are those obtained by 
refining T. For any ASC p and any 1 < i < 4 we denote 

Vi(p) := {T E Vt] p(T) does not hold}. (24) 



The first point of Definition 2.4 implies that Vi(p) is a (possibly empty) subtree of Vf. any 



T' E Vi{T) is either the root Tj, or the child of another T E 7^. The second point of the same 



definition, combined with (22), implies that Vi(p) is finite 



The finiteness of the trees Vi(p), implies that the refinement procedure ends after a finite 



number of steps. As already observed right after Definition 2.2 the resulting collection of triangles 
is automatically a mesh satisfying Points (i) and (ii) of Definition 

As observed in the second point of Lemma |2,3[ the collection of elementary triangles is the 
disjoint union of the trees ("Pj)i<j<4. Thus the subtrees Vi{p) form a partition of the set £{p). 
Recalling that the number of leaves of a binary tree is one plus the number of its inner nodes, 
we obtain 

#(T(p)) = £ (1 + #(7>i(p))) = 4 + #(£ (p)), (25) 

l<i<4 

which concludes the proof of the first point. 

The implication p =>• p' of two ASCs is equivalent to the reverse implication of the negations 
->p <= -ip', and thus to the inclusion £{p) D £{p')- If p => p' we thus obtain #(£(p)) > #(<5(p')), 
and therefore #(7~(p)) > #(T(p')), which establishes Point 2. 

For any two ASCs p,p', we have £(pAp') = £{p) U£(p'). Hence #(£(p Ap')) < #{£(p)) + 
#(<5(p')), and therefore #(T(p A p')) < #(T(p)) + #(T(p')) — 4, which concludes the proof of 
this proposition. □ 

We establish in the following lemma a first non-trivial estimate of the mesh cardinality 
#(T(p)), in terms of the constant s p associated to the ASC p. 

Lemma 2.7. • For any s E [l,oo[ one has 

£ ^< 8(1 + Ins). (26) 
0<||«|[<s 

• For each u E 72? denote 

£u := i v e ll n ll < IMI> < < s p , det(u,v) = 1}, 

and define £~ likewise, to the exception of the last constraint which is replaced with 
det(u,v) = -1. Then < s p /||n|| 2 , fore E {+,-}. 



14 



There exists a constant C such that for any ASC p, with associated constant s p > 1, one 
has 

#(T(p))<Cs p (l + lns p ). (27) 



Proof. We first establish (26), and for that purpose we introduce the sup-norm || • ||oo on H 2 
defined by [|(a;, 2/)[|oo := max{|x|, |y|}. Clearly || 

^||oo — II ^|| fo^ ^ ^ 1R 2 . For each k G 
there exists precisely (2k + l) 2 elements u G H? such that ||m||oo < k. Hence for each integer 
k > 1 there exists precisely (2k + l) 2 — (2k — l) 2 = 8k elements u G 2Z 2 such that ||ti||oo = k. 
Therefore 

E isp* E ht- E pS8d+M, 

«ezz 2 11 11 «ea 2 11 "°° o<fc< s 

0<||«||<s 0<||m|| oo <s 



which concludes the proof of ( 26 ) 



We next turn to the proof of the second point, and for that purpose we consider a fixed 
u G 7L 2 such that is non-empty. Let v G £+ be such that the scalar product (it, is 
minimal. For any v' G one has det(u, v' — v) = 1 — 1 = 0, hence u ' = v + An for some A G JR. 
Since u,v,v' G 2Z 2 , and since « has coprime coordinates (recall that det(u,v) = 1), the scalar A 
must be an integer. We thus have 

IM| 2 < (u, v) < (u, v 1 ) = (u, v) + A||u|| 2 < s p , 



where we used the first point of Lemma 2.3 for the leftmost inequality. Hence < A < s p /||u|| — 
1, and therefore #(<f^) < s p /||u|| 2 . Estimating £~ likewise, we conclude the proof of the second 
point. 

Identifying an elementary triangle to its pair (u, v) of non-zero vertices, ordered by increasing 
norm, we find that £(p) \ 7o injects in the union of the £^, for u G 2Z 2 and e G {+,—}. 
Furthermore £^ is empty for ||u|| > y/sp~, since using the second point of the proposition we find 
that its cardinal is strictly less than one. Hence using the first point of this proposition 

#(£(p)) <#(%)+ £ (#(£+) + #(0)< 4 + E |^<4 + 16 Sp (l + ln Sp ). 

ueK 2 0<\\u\\<y/^p 



Recalling that #(T(p)) = 4 + #(£ (p)), we conclude the proof of this proposition. □ 
2.2 Mesh associated to an asymmetric norm 



We describe in Proposition 2.9 of this section construction of a F-acute mesh T(F) for each 
asymmetric norm F. 

Our first lemma introduces a tool that will be frequently used in the subsequent analysis: 
the approximation of an arbitrary asymmetric norm by smooth ones. For each G 1R we denote 

e e := (cos 0, sin 0). (28) 

Lemma 2.8. For any asymmetric norm F on M 2 , there exists a sequence (F n ) n >i of asymmetric 
norms such that 

• F n —7- F locally uniformly as n —> oo. 

• F n G C°°(i?? 2 \{0}). 

• n(F n ) < k(F) for alln>\. 
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Figure 6: Unit ball of an asymmetric norm F of anisotropy ratio k(F) = 20 (top left). 
Generation of T(pf) by recursive bisection (bottom, left to right). The non-zero vertices of 
colored triangles do not form an F-acute angle, hence these triangles are refined. Colored 



triangles constitute the set £(pp), see definition 2.5 they form the inner nodes of four finite 



binary trees of triangles, while the elements of T(pf) are the leaves. 

• If F is symmetric, then so is F n for all n > 1. 

Proof. We define the asymmetric norm F n through polar coordinates and by convolution; for 
each r > and each ip 6 H, 



F n {re v ) := r [ F{e e )n n {v ~ B)M. 



'JR. 

We denoted by fi n the mollifier fi n (9) := n/j,(nO), for all n > 1, where fj,(0) := e _£,2 / v / 7r, for all 
9 G JR. The four announced properties are immediate. □ 

We presented in the introduction of this paper the construction of a mesh T(F), associated 
to each asymmetric norm F. This definition is tied to the mesh generation method by recursive 
refinement presented in the previous subsection, since we claim that 

T(F) = T( PF ), 

where for an elementary triangle T of non-zero vertices u,v, the predicate value Pf(T) stands 
for the test u u,v form a F-acute angle", see the next proposition. Indeed, denote by (7fc)^L 
the elementary triangles defined by the consecutive pairs u, v of vectors subject to the test "If 
u,v form a -F-acute angle", in the construction of T(F). These triangles, and their indices, are 
shown on Figure [6j The sequence (Tk)^ =0 constitutes an in-order transversal of the four binary 
trees in £{pf) U T(pf), see again Figure [6] and the proof of Proposition 2.6 



Let us observe that Point (i) and (ii) of Definition [lTJ of F-acute meshes, hold by construe 



tion for any mesh of the form T(p), where p is an ASC, as observed right above Lemma 2.3 On 
the other hand the predicate pf is designed so as to enforce Point (iii) of this definition. 
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Proposition 2.9. Let F be an asymmetric norm on M, . 

• Two vectors u,»£ M 2 \ {0} form an F -acute angle if 

k{F)svd.\<{u,v)\ < 1. 

• The predicate pf defined for any elementary triangle T by 

Pf(T) holds if and only if the non-zero vertices of T form a F -acute angle, 



(29) 



(30) 



is an ASC, with associated constant s PF < \J k(F) 2 — 1. 
Any vertex u ofT(F) satisfies \\u\\ < 2n{F). 



Proof. In order to establish (29), we restrict in a first time our attention to asymmetric norms 
F which are smooth: F G C^(R 2 \ {0}). In that case Proposition |3,6| (below, but proved 



independently) shows in (J42J) that k(F) cos <{cq, VF(eg)) < 1, for any 9 € K. Since F is 
homogeneous, one has V-F(w) = VF(Xu) for any A > and any it 6 H 2 \ {0}, and therefore 

k(F) cos <(u,VF(«)) < 1. 



Assuming (29) we thus obtain 



\<(v,VF(u))\ < \<(v,u)\ + \<(u,VF(u))\ < arcsin(l/K(F)) + arccos(l/K(F)) = vr/2, 



and therefore (v,X7F(u)) > 0. Likewise (u, ~VF(v)) > 0, hence using Point 1 of Lemma 1.2 we 
conclude that the vectors u, v form an F-acute angle. 

Now let us consider an arbitrary, possibly non-smooth, asymmetric norm F, two vectors u, v 
satisfying (29), and a sequence (-F n ) n >o of asymmetric norms as described in Lemma 2.8 It 



follows from the above argument that the vectors u, v form an i^-acute angle for each n > 0. 



Since Definition 1.1 of F-acuteness only involves non-strict inequalities, we obtain taking the 
limit that u, v form an F-acute angle. 

We next turn to the proof of the second point of this proposition. Consider two vectors 
u, v 6 H 2 \ {0} which form an i^-acute angle. For each 5 > we obtain 

F(u + v + 5u) = F((l + 5)(u + v)-5v) 

> (l + 8)F(u + v)-8F(v) 

= F(u + v) + 5(F(u + v)-F(v)) 

> F(u + v), 

where we used the triangular inequality in the second line, and the fact that u and v form an 
F'-acute angle in the last. For the same reason F(u + 5(u + v)) = (1 + 5)F(u + 5v/(l + 8)) > 
(1 + 5)F(u) > F(u). Therefore u and u + v form an F-acute angle, and likewise v and u + v 
form an acute angle. As a result, if the predicate pf holds for an elementary triangle, then 
it also holds for its children. This establishes the heredity property in Definition |2.4| For the 
finiteness property we consider an elementary triangle T, of non-zero vertices u, v, such that 
s(T) > ^Jk(F) 2 - 1. It follows from (21 ) that sin \<(u, v)\ < 1/k(F), hence the first part of this 



proposition shows that u,v form an F-acute angle. Thus Pf{T) holds, which establishes the 



finiteness property of Definition 2.4, and concludes the proof of the second point. 
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We turn to the third point of this proposition. A triangle T £ T(F) either belongs to To, or 
is the child of a triangle T" in £(pf) which does not satisfy the predicate pp. If the former case 
there is nothing to prove, while in the latter case the non-zero vertices u, v of T" satisfy 

max{||u|| 2 , \\v\\ 2 } < \\u\\ 2 \\v \\ 2 = s(T) 2 + 1 < s 2 pF + 1 < k(F) 2 . 



We used the fact that min{||u||, \\v\\} > 1, since these vectors have integer coordinates, and (21). 
Thus \\u\\ and ||u|| are bounded by k(F), and therefore \\u + v\\ < 2k{F). This concludes the 
proof since the non-zero vertices of the triangle T belong to {u, u + v, v}. □ 



At this point, we can establish the worst case analysis presented in Proposition 1.5 except 



for (17) which is proved later in Corollary 3.10 



Corollary 2.10. • Let k > 1, let p K be the predicate "s(T) > k", and let T K := T(p K ). 
Then T K is F -acute for any asymmetric norm F such that k(F) < k, and furthermore 
#(T(F)) <#(T K ) <C7c(l + bi/c). 

• For each r > 1 let F T be the anisotropic euclidean norm defined by the positive definite 

matrix M T := ( * ^ T 2 \ Then \ K ( F r) - 2t| < 1 and #(T(F T )) > 6 + 2[rJ . 

Proof. We begin with the first point. For any asymmetric norm F such that k(F) < k, we 
have s PF < \J k(F) 2 — 1 < \/ k 2 — 1 < k. Hence p K =4> pf, which implies simultaneously that 
#(T(-F)) < #(7^) (using Point 2 of Proposition 2.6) and that T K is .F-acute (since pf holds for 



all the elements of T K ). The estimate on #(7^) was proved in Lemma 2.7 



We turn to the proof of the second point. The 2x2 symmetric matrix M T is positive definite 
since its trace and determinant are both positive. Denoting by < A 2 < fi 2 the eigenvalues of 
M T , where A and /i are positive, one has 

K (F T ) = ^, Tr(M r ) = A 2 + ii 2 = 2t 2 + 1, det(M T ) = AV = 2r 2 - r 2 = r 2 . 
A 

Hence denoting k := &{F T ) 

l +ft= A + ^ = ^ L = 2^ ±i = 2T+ l 
k /i A vdetM T r r 

Therefore \k — 2r| = |r _1 — < 1, since k > 1 and r > 1. 

We next observe that the elementary triangle of vertices (1, 0) and (— r, 1) is not i^V-acute for 
< r < r, since ((r, — 1), M r (l, 0)) = ((r, —1), (1, r)) = r — r < 0. Considering these triangles 
and the symmetric ones with respect to the origin, we obtain 2(1+ [tJ) non F T - acute elementary 
triangles, hence #(T(F T )) = 4 + £(pp T ) > 6 + 2[tJ using (23), which concludes the proof of the 
first point. □ 

In the next section, we estimate the cardinality of T(F), for asymmetric norms F which are 
smooth on H 2 \ {0}. These results are transferred to arbitrary asymmetric norms, using the 
approximation result Lemma 2.8 and the following lemma. 

Lemma 2.11. Let F be an asymmetric norm, and let (F n ) n >o be a sequence of asymmetric 
norms such that F n — >■ F locally uniformly as n — ^ oo. Then 

#(T(F)) < liminf#(T(F n )), (31) 

n— >oo 

r#(T(F e ))de < hminf r#(T(F*))d9. (32) 
Jo n ^°° Jo 
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Proof. To avoid notational clutter, we denote £(F) := £{pp), and £{F n ) := £{Pp n ) 
inition 2.5 If an elementary triangle T belongs to £(F), then it belongs to £(F n ) 
sufficiently large, since .F-acuteness is a closed condition, see Definition 1.1 Hence 



see Def- 
for all n 



£{f) c u n *m 



N>0n>N 

This immediately implies that #(£(F)) < limmfn-^oo #(£(F n )), by applying Fatou's lemma to 
the characteristic functions of £ (F) and £(F n ). Inequality ( |31| ) then follows from the identity 
#(T(J0)=4 + #(f(F)),seeg3j. 

The second estimate (32) immediately follows from the first one (31 ), by observing F® — > F e 
locally uniformly as n — > oo for any 8 G [0, 2ir], and applying Fatou's lemma on this interval. □ 



3 Average Complexity 

This section is devoted to the estimate of the size #(T(F)) of the stencils using in the FM-ASR, 
or more precisely of the average value of #(7~(F 8 )), as 9 G [0,2-7r]. Estimates are obtained for 
increasingly general types of (asymmetric) norms F. We begin with anisotropic euclidean norms 
in the first subsection, then turn to symmetric norms in the second, and finish with asymmetric 
norms in the third. Each subsection builds on the estimate of the former one, hence they are 
not independent. 



3.1 Anisotropic euclidean norms 

Our first lemma shows that the triangles refined during the construction of T{F), where F is an 
anisotropic euclidean norm defined by a symmetric positive definite matrix M, are all aligned 
with the eigenspace associated to the small eigenvalue of M, see also Figure [7} 

Lemma 3.1. • Let F be an anisotropic euclidean norm, given by a matrix M G S 1 ^. If 
the non-zero vertices of an elementary triangle T do not form an F -acute angle, then T 
contains a (non-zero) eigenvector for the smallest eigenvalue of M . 

• For any anisotropic euclidean norm F one has ff{T{F)) < 6 + 2k{F). 

Proof. In order to establish the first point, we denote by < A < \i the eigenvalues of M, and by 
e a normalized eigenvector of M associated to the eigenvalue A. We denote by u, v the non-zero 
vertices of T. We have 

> {u, Mv) = X(u, e)(v, e) + /zdet(it, e) det(f , e) = X(u, v) + (fi— A) det(u, e) det(u, e), 

where we used the identity (u, v) = (u, e){v, e)+det(u, e) det(v, e). Since u, v do not form an acute 



angle, we have (u,Mv) < using Point 2 of Lemma 1.2 On the other hand (u, v) > 0, hence 
det(u, e) det(t>, e) < 0, and therefore det(u, e) and det(t> , e) are non-zero and have opposite signs. 
It follows by continuity (or linearity) that there exists t G (0, 1) such that det(tu+(l— t)v, e) = 0, 
which concludes the proof of this point. 

We next turn to the proof of second point, and for that purpose we adopt the notations 



of Proposition 2.6 and consider the four trees (Vi(pF))i<i<A- It follows from the first point of 
this lemma that two of these trees are empty, and that the other two have a single branch, see 
also F igure |7| The number of elements of these single branched trees is bounded by 1 + sp F 



1 + ^ k(F) 2 - 1 < 1 + k(F), using ft22n and the second point of Proposition 2.9l We finally 



obtain using (25) that #(T(F)) <4 + 2x0 + 2(l + k(F)) = 6 + 2k{F) which concludes the 



proof. □ 
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Figure 7: Unit ball {u; F(u) < 1} of a norm F given by a positive definite matrix M, of 
anisotropy ratio k(F) = 8 (left). Eigenspace associated to the small eigenvalue of M (dotted 
line). Generation of T(F) by recursive bisection (second left to right). All refined triangles 
(colored) contain an eigenvector associated to the small eigenvalue of M in their interior. 



Definition 3.2 (The following definitions are restricted to this section). We consider a fixed 
constant k > 1, and denote by F the anisotropic euclidean norm defined by the diagonal matrix 
D of entries (k~ 1 , k), in such way that k(F) = k. For each 9 G [0, 2tt] the norm F e is also of 
anisotropic euclidean type, and defined by the matrix Mg = RgDRj . 

• We denote by £g, 9 G [0,2n], the collection of elementary triangles which non-zero vertices 
do not form a Fg-acute angle. In other words Eg := £(p F e). 

• For each elementary triangle T, we define It '■= {9 £ [0, 2ir]; T S £g}. 



It follows from (23), that for any 9 G [0, 2tt] one has 

#(T(F 9 )) = 4 + #(£ e ). (33) 
Furthermore, we have by construction 

#(£ e )d9 = T\I T \, (34) 
Jo T 

where T ranges over all elementary triangles, and \It\ denotes the Lebesgue measure of It- 

Proposition 3.3. For any fixed u G 2Z 2 \ {0}, we denote by A u the collection of all elementary 
triangles T , containing u as a vertex, and such that the other non-zero vertex v satisfies \\u\\ < 
\\v\\. We have 

E i*i <- (^) 

and furthermore this sum equals if \\u\\ > ^[k. 

Proof. Let T be an elementary triangle such that It ^ 0, and let u, v be the non-zero vertices 
of T, with ||u|| < \\v\\. Since n(F e ) = k for any 9 G 1R, we obtain using (29) and (20) 

- < sm\<l(u,v)\ = „ „ < t-^, (36) 
K \\u\\\\v\\ nr 
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and therefore ||u|| < \fn. It follows as announced that (35) is zero if > 

It follows from Lemma (3.1) that tu+ (1 — t)v is proportional to eg, for some t G]0, 1[. Hence 

\I T \ < 2|<(«,«)| = 2arcsin ( „ „ ) < „ ^ „ , (37) 

V \\ u \\ \\ V \\J \\ u \\\\ v \\ 

where we used the concavity estimate sin(7rx/2) > x for x G [0, 1]. 

We denote by v e , for e G {—1, 1}, the element of 7L 2 for which the scalar product (u,v £ ) is 
non-negative and minimal, under the constraint that det(u,v £ ) = e and \\v E \\ > \\u\\. 

If det(u,v) = e, then v — v £ = Xu for some A G R. Observing that u,v,v £ have integer 
coordinates, and that u has coprime coordinates, since det(it, v) = 1, we obtain that A is an 
integer. The scalar A is non-negative since (u,v £ ) < (u,v) = (u,v £ ) + A||n|| 2 . Last we observe 
using (36) that k > \\u\\ \\v\\ > (u, v) = (u, v £ + An) > A||u|| 2 , hence A < k/||u|| 2 < k. 

We have \\v £ \\ > \\u\\ by construction, and \\v £ + Xu\\ > A||u|| for any A > 1, since (u,v £ ) > 0. 
Therefore, recalling ( [37] ), 

7r 2n < 27r(2 + lnK) 



£i'H< E E 



Jitll llug + A«|| ~^ llnll 2 max{A, 1) ||n|| 2 

T&A U ££ {i_i} 0<A<ft 6 11 0<k<K 11 11 11 11 

which concludes the proof of this proposition. □ 



We next establish Theorem 1.4, in the special case of anisotropic euclidean norms. 

Corollary 3.4. There exists a constant C , such that for any anisotropic euclidean norm G one 
has 

,-2-k 



r ~#(T(G e ))<C(l + \n K (G)) 2 . 
Jo 



Proof. It is sufficient to prove this result for the specific norm F introduced in Definition 3.2 



since any anisotropic euclidean norm has this form, up to a rotation and a multiplication by a 
positive scalar. 

The sum (34) of the interval lengths \It\ associated to all elementary triangles T, can be 
bounded as follows 

El J H< E E |/r|<C(l + ln«) E li <8C(l + ln K ) 2 , 

T ue7L 2 T&A u «eS 2 \{0} " " 

l«ll<\A: 



where we used (26) for the last inequality. Combining this estimate with (33) we conclude the 
proof 

r>27r 

□ 



/ #(T(F e )) <4 + 8C(l + hiK) 2 . 
J o 



3.2 Symmetric norms 



In this section and the following one, we denote by 5 the collection of asymmetric norms which 
are continuously differentiable outside of the origin. To each F G $ we attach a 27r-periodic map 
ipp : H —7-] — 7r/2,7r/2[, introduced in the following definition, which describes the direction of 
its gradient. See Figure [8] (left) for an illustration and (center) for two examples. 



Definition 3.5. For any F G "$ and any 9 G M we define the angle (20) 

Vf{0) :=<(e e ,VF(e e )). 
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Figure 8: Illustration of Definition 3.5 (left), ^>f{9) < in this example. Graph of i -Pf{9) (center) 
for an anisotropic euclidean norm given by a diagona l ma trix (plain) , and the asymmetric norm 
■sj x 1 + y 2 — 0.9x (dashed). Notations of Proposition 3.6 



Since the asymmetric norm F is 1-homogeneous, we have {eg,X?F(eg)) = F{eg) > for all 
9 G H, hence 

p F (0)€]-?r/2,7r/2[. (38) 
The composition of F with a rotation, corresponds to the composition of ipF with a translation: 



ip F e = ip F (- ~ 0). 



(39) 



Note that (pp is 7r-periodic if F is symmetric, and odd if F(x, y) = F(x, —y) for all x, y G H. 
The next proposition establishes the two most noticeable properties of ifF aside from its 



periodicity: its integrals are bounded (40), and it obeys a semi-Lipschitz regularity property 



(41). 



Proposition 3.6. For any F G 5 and any 9 G M, one has ^lnF(eg) = t&n(pp(9). Hence for 
any h > 



tanyj^ 



< lnK(F). 



Furthermore ipp is right-Lipschitz, and \(pp\ is bounded strictly away from ir/2: 



(p F (0 + h) 
cos <pf(0) 



> 
> 



Vf(0) - h, 
1/k(F). 



(40) 



(41) 
(42) 



Proof We define r G C 1 (R,1R+) by r(6>) := l/F(e e ), for all 6> G M. This quantity is illustrated 
on Figure [8] (right), as well as the vectors ejf and Vi ? (eg)" L , where u 1 - denotes the rotation by 
7r/2 of a vector u G H 2 . We immediately obtain the Taylor development 

r(9 + h) = r(6) + hr(6) tsn.(-tp F (0)) + o(h), 

for any fixed fl £ R, and for small h. In other words r'(#) = — r(9) tainpp(6), equivalently 
InF(ee) = tan(pp(9). The left hand side of (40) therefore equals \lnF(eg) — lnF(t 



de 



[ee+h) 



which as announced is bounded by Itik(F), by definition of the anisotropy ratio Q. 

Let i/j(0) := 6 + (fp(0), for all 9 G JR. By construction VF(eg) is positively proportional to 



for all 



G H. The vectors e^^ and are respectively the unit normal and the unit 



tangent to the set Bp := {z G M 2 ; F(z) < 1}, in the direction eg. Since Bp is convex, the 
derivative of the tangent vector ^g-e i p(9) 1 ~ = — V''(#) e ?/j(6>) is negatively proportional to normal 
This shows that ip'{0) > 0, for all 9 G M, hence that ifi is non-decreasing. Recalling that 
i -Pf{9) = ip{0) — 9, we conclude that cpp is the difference of a non-decreasing function and a 



1-Lipschitz function, which establishes (41). 
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For the last inequality, we fix 6 and first assume that p := Pf(6) > 0. We obtain using (40) 

rO+tp rp 



ro+p rp 
In k(F) > / tan p p > / tan 
Jo Jo 



(p — u)du = — ln(cos tp), 



hence cos(p) > l/n(F), as announced. If p < 0, a similar argument involving the integral on 
[0 + p, 6] yields the same estimate, which concludes the proof of this proposition. □ 

We rephrase in the next lemma a geometrical property, on the gradients of a family of 
asymmetric norms, into inequalities between the attached functions. 

Lemma 3.7. Let F, F±, ■ ■ ■ , F r € and let u £ M 2 . The following are equivalent: 

• There exists an, ■ ■ ■ ,a r £ 1R+ such that 

VF(u) = <*iVFi(u)- (43) 

l<j<r 

• Let 9 G M be such that u and eg are positively collinear. Then 

min{ m (0),-- - , Vf M}<M0)<^{<PfM>"- ,<PF r (P)} (44) 

Proof. Since F is 1-homogeneous, we have VF(u) = X7F(u/\\u\\) = VF(eg) (likewise V-Fj(ii) = 
V-Fi(eg)). We denote v := AVF(u) (resp. Vi := AjVFj(u)), where the positive scalar A (resp. 
Aj) is chosen so that (u,v) = 1 (resp. (u,Vi) = 1). We introduce the angles p := Pf(0) 



(resp. pi := pp^O)), which belong to ] — 7r/2,7r/2[, see (38), and we observe that tan</? 
det(u,v)/ (u,v) = det(u,v) (resp. tan^j = det(w, Vi)). 

Assuming (43), and denoting := (A/Aj)aj > 0, we have v = Yli=i fii v i an d therefore 

1 = (U,V)= ^ Pi{u,Vi) = ^ Pi, 
l<i<r l<i<r 

tan p = det(u, v) = det(u, V{) = ftitanpi. 



KKr Ki<r 



This shows that t&np is a weighted average of the (tan Pi) r i= i, which implies (44) since tan is 
increasing on ] — 7r/2,7r/2[. 



Conversely, if (44) holds, we may assume without loss of generality that pi < p < p2- We 
thus have tany?i < tan p < tan p2, hence there exists barycentric coefficients /3i,@2 £ H+, 
Pi + P2 = 1, such that tan p = Pi tan^x + /3 2 tan p2- Setting /3 3 = • • • = /3 r = 0, we thus have 
tan 93 = ^[ =1 /3j tan 9?j. Defining V = ^ ftvj, we obtain proceeding as above (u, v) = 1 = 
{u, V) and det(u, v) = tanp = det(u,V), hence v = V. Denoting on := (Aj/A)/3j we obtain 



VF(u) = Ya=i «tVi4(ii), which establishes (43) and concludes the proof. □ 



We emphasize the next proposition, which is a central component of our strategy. Consider a 
"complex" asymmetric norm F, and "simpler" norms {Fi) r i=1 , say of anisotropic euclidean type. 



Assume that ipp is bounded in the sense of (45) by the p Fi . The following result shows that 



#{T{F)) can be estimated in terms of the #(T(-Fj)), for which efficient bounds were developped 
in the previous subsection. 
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Proposition 3.8. Let F,F±, - ■ ■ ,F r 6 $ be such that everywhere on M 



mm{<p Fl , ■ ■ ■ , <PF r } < < raax{ip Fl , • • • , ip Fr }. 



(45) 



Then 



#(T(F))< #(nm, <™d r#(T(F 9 ))d0< ]T r#(T(F t 

Kr<r Ki<r 



))d6. (46) 



Proof. We begin with the proof of an intermediate result: if (45) holds, then we have the 
implication of ASCs 

PFi A • • • A pF r Pf ■ (47) 

Indeed let T be an elementary triangle, of non-zero vertices u, v. Let 9 U ,9 V G 1R be such that 
eo u ,eg v are respectively positively proportional to u, v. Assume that {pp x A • • • f\pF r )(T) holds, 
which means that u, v form an Fj-acute angle for all 1 < i < r. Using Lemma 



L2 
3~7 



we find that 
we find that 



{u, V-Fj(f)) > and (v,VFi(u)) > 0, for all 1 < % < r. Using (45) and Lemma 
V-F(ii) is a linear sum with non-negative coefficients of V-Fj(u), 1 < i < r, hence {v , X7F(u)) > 0. 
Likewise {u,VF(v)) > 0. Using again Lemma 1.2 we obtain that u, v form a -F-acute angle, 



hence Pf{T) holds. This concludes the proof of (47) 



The left part of (46) immediately follows from (47) and Point 3 of Proposition 2.6 Due to the 
translation invariance ( 39 ) , ( 45 ) is equivalent to minj^e , • • • , tp F e } < (p F e < maxj^e , • • • , (f F e } 
for any 9 G R, and thus implies ff{T{F e )) < Yli<i< r ))• Integrating over [0, 2-k] we ob- 

tain the right part of (46), which concludes the proof. □ 



The next technical lemma investigates the periodic function attached to an anisotropic eu- 
clidean norm. This is a prerequisite if one wants to construct a well chosen family (i ? j)[ =1 of 



such norms which satisfies (45), given an asymmetric norm F of interest. 



Lemma 3.9. • The anisotropic euclidean norm F defined by the diagonal matrix of entries 
(k ,k), where k > 1, satisfies k(F) = k. The function ip F attains its maximum at 
9 K := arctan(K _1 ) ; which is 



1 Pf{9 k ) = arctan [(7 



/2] 



• There exists a finite number anisotropic euclidean norms F±,--- ,F r , such that on M 

mm{(f Fl ,- ■ ■ ,(f Fr } < -tt/4 and ir/4 < max{^ Fl , • • • , fF r }- 

Proof. Let D be the diagonal matrix of entries (1/k, k). We have F(u) 2 = (u, Du) for any u £ 
1R 2 , hence F(u)VF(u) = Du for any u ^ 0, by differentiation. It follows that F(eg)VF(eg) = 
(k _1 cos 9, Ksin#), for each PgR, and therefore 



tamp F {9) 



K — 1/k 



det(ee, VF{ee)) (k — 1/k) cos#sin6> 
(eg,VF{ee)) n^ 1 cos 2 9 + Ksin 2 9 (Ktan6*) _1 + Ktan6> : 



where the right hand side equals by convention if 9 is a multiple of tt/2. The maximum value 
of this right hand side is attained when its denominator is positive and minimal, that is when 
«tan# = 1. Thus the maximum value of t&n(p F is (k — 1/k)/2, attained at 9 K , which concludes 
the proof of the first point of this proposition. 
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We next turn to the proof of the second point, and for that purpose we choose k sufficiently 
large (k = 13 is fine) in such way that arctan [wc — k^ 1 ) /2] > tt/4 + tt/5. We denote by F the 
norm associated to the diagonal matrix of entries (1/k,k). 

Since the norm F is symmetric, the function ip F is 7r-periodic. Since F is defined by a diagonal 
matrix we have F(x, y) = F(x, —y) for all x, y G M, and therefore -F is odd. Furthermore for all 



& G + tt/ 5 ] we have ^(9) > <£ F {6 K ) - tt/5 > tt/4, using fl4l] ). Finally for all n € 2Z 

> tt/4 on [n7r + # re , nir + # K + vr/5], and ipp < — 7r/4 on [mr — K — tt/5, n7r — 
We choose r := 5, and introduce the anisotropic euclidean norms i 7 ^ := F k7T ^ 5 , for 1 < fc < 5. 



Using the translation invariance (39 ) we see that the sets on which (p Fj > tt/4 (resp. tp Fi < —tt/4), 



1 < j < r, contain intervals which cover the whole line M. This concludes the proof. □ 

For any asymmetric norm F and any A £ GL2, we denote hy F o A the asymmetric norm 
defined by 

F oA(u) := F(Au). 

Clearly FoA is symmetric (resp. is an element of J 7 , resp. is of anisotropic euclidean type) if and 
only is that is the case of F. We denote k(A) := ||^4|| ||^4 _1 ||, and point out that || Au||/|| Au|| < 
k(A) for any u, v 6 M 2 such that ||w|| = \\v\\ = 1. It easily follows that 

max {«pj'^} - k(FoA) - < F ^< A y ( 4§ ) 

Choosing A = Rj , for some £ R, we recover in particular that n(F e ) = k(F), for all 9 €zR. 
We establish in the next corollary the average case and the worst case estimates for ^(F(F)), 



where F is an arbitrary symmetric norm, which were announced in Theorem 1.4 and Proposition 



1.5 respectively. 



Corollary 3.10. • For each symmetric norm G there exists A G GL2 such that F := Go A 
satisfies k(F) < y/2. If G G J, we therefore have on M 

- tt/4: <(f F < tt/4 (49) 
• There exists a constant C such that for any symmetric norm G one has 

#(T(G)) <Ck(G), and / #(T(G e )) d9 < C(l + In k(G)) 2 . (50) 

Jo 

Proof. The minimum value of k(G o A), for ^4 G GL2, is attained and coincides with the multi- 
plicative Banach-Mazur distance d(G, \\ ■ ||) between the norm G and the euclidean norm || • ||. 
It follows from John's theorem on the maximal ellipsoid that this distance is bounded by y/2 for 
2-dimensional norms {\fd for d-dimensional norms). Hence we may choose A G GL2 such that 
F := G o A satisfies k(F) < If G € we thus have cos(<pp(9)) < I/k(F) < l/y/2, for all 



6 G IR, using Proposition 3.6 Finally |</3p(0)| < arccos(l/v2) = tt/4, which concludes the proof 
of the first point. 

We keep the same notations, and turn to the proof of the second point. In view of Lemmas 



2.8 



3.9 



and |2.li| we can assume that G G 5- Let F\, ■ ■ ■ , F r be as in the second point of Lemma 
and let Gi := F{ o A -1 , for 1 < i < r. We have by construction 

mm{ip Fl ,--- ,ip Fr } < -tt/4 < ip F < tt/4 < max{(/? Fl , • • • ,(p Fr }, 
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on 1R. Hence for each u G 1R 2 \ {0} there exists, using Lemma 3.7, non-negative coefficients 
£*!,••■ ,a r such that: 



A T VG(Au) = VF(u) = Y aiVFi(u) = A T I Y ^G^Au) 



(51) 



It follows that VG(n) is a linear combination with non-negative coefficients of the VGj(t>), 

we conclude 



1 < i < r, for any v G H 2 \ {0} (choose u = A 1 v in (51 )). Using again Lemma 



3.7 



that min{ip Gl , ■ ■ ■ , <p Gr } < <p G < max{(p Gl , ■ ■ ■ , fc r } ° n R. 
Hence 



#(T(GQ) < Y #V(Gi)) <C Y K (Gi), (52) 

l<i<r l<i<r 

#(T(G 8 )) d9 < Y, r rt^G*)) < C Y (1 + ln K (^)) 2 , (53) 

Ki<r l<j<r 



where we used Proposition 3.8 for the first inequality, of both lines, and for the second inequality 



Lemma 3.1 in the first line, and Corollary 3.4 in the second line. 
On the other hand we obtain using (48), with k(A) := ||.A||||.A _1 | 



k(G,) = o A) < K (Fi) K {A) and «(G) > = ^ > 



thus k(Gj) < V2K(Fi)K(G). Combining this inequality with (52) and (53), we conclude the 
proof of d50b. □ 



3.3 Asymmetric norms 

We estimate in this section the average cardinality of #(T(F S )), € [0, 27r], for an arbitrary 
asymmetric norm F on H 2 . Our strategy is similar to the case of symmetric norms, presented in 
the previous subsection: we construct a family F\, ■ ■ ■ ,F r of anisotropic euclidean norms such 



that mm{(fF 1 , • • • , <fF r } < <-Pf < max{(pp 1 , • • • , <pF r }, and we use Proposition 3. 



The construction the {Fi) r i=l is however more subtle than in the previous section, and the 
integer r grows logarithmically with n(F). The following technical lemma, illustrated on Figure 
[9j is our first step. 

Lemma 3.11. Let < e < vr/6, and let $ G C°([0, 2vr], ]0, vr[) be such that (i) $ > 2e, (ii) 
<3? — Id is non-increasing, and (Hi) for all < t < t' < 2tt one has 



ft' 

/ cotan<&(s) ds < | In e\. 



(54) 



Let (t n )o< n <N be the finite sequence of elements o/[0, 27r] recursively defined as follows: to := 0, 
and i n +i is the largest solution in [0, 2tt] of 



<f>(t) = e + (t-t n ) 

if one exists. Otherwise the sequence ends. 

Then N < C\ lne|, for some absolue constant C (independent of e and 



(55) 
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Figure 9: Illustration of Lemma 3.11 Curve : <£. Broken line : construction of to,-- - ,ij\r- 



The intersections of the broken line with the curve correspond to the equations (55). There is 
one vertical segment at each abscis sa t n , < n < N. The chosen $ is ir/2 

3.121 with F(x,y) 



consistently with Proposition 



V* 2 



ipp, or 7r/2 — tfp* 



y 



2 - 0.98x. 



Proof. Since t \-> &(t) — t is continuous and non-increasing, the solutions of (55) are either the 
empty set, a singleton, or a closed interval. We denote (ft n := <&(t n ) for all < n < N, and 
observe that <j) n > 2e by hypothesis (ii). We also define for < n < N 

<5n := tn+l — t n = 4>n+l ~ £>2e — E = £. 

In particular t n > tie for < n < N, and therefore N < 2tt/e: the sequence is finite as 
announced. We need however to estimate the integer iV more carefully. The proof below 
heuristically distinguishes for each segment [t n , i„,+i] three possible cases. The segment is either 
(i) long, if t n+ \ — t n > 7r/6, or (ii) such that 4> n +i <fi n , or (in) such that cotan $ is 

significant. A bound on the number of segments of type (i) is first obtained, and then a bound 
on the number of consecutive segments of type (ii) or (iii) between to segments of type (i) . 
We introduce two collections of integers 

E := {2, • • • , JV - 1}, E + := {n G E; 5 n _i + 5 n < tt/6}. 

We first observe that 

(it/6)#(E\E + ) < (Sn-l + 5 n )<2 $n = 2(t N -t ) <4tt, 

n&E\E + 0<n<N 

and therefore #(E \ E + ) < cq ■= 24. We next estimate the cardinality of E + . Consider an 
arbitrary n G E+, we obtain recalling that $ — Id is non- increasing 

[tn+i r$n /sm((h +5 )\ 

/ cotan dt > / cotan(</) n + 1) dt = In — — — . 

Jt n JO V Sm ^n ) 

We next define and estimate a quantity e n , attached to each n G E + 

«, :_ ,„ t^L\ + 2 /•"■« cotan , > ln ( ) . (56) 

Vsin(0 n+ i)y J tn \sm((t> n ) sm(0 n+ i) J 

We have sin(^> n ) < (f> n , and sin(0 n+ i) < n +i = e + S n < 2S n , since 5 n ^ e. We recall on the 
other hand that sin(x) > px for all x G [0, 7r/3], by concavity of the sine function on this interval, 
with p := sin(7r/3)/(7r/3). Observing that <j) n + 5 n = £ + <5 n _i + S n < ir/Q + tt/Q = tt/3, we thus 
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obtain sin((^ ri -|- dn) ^ p((fr n -\- S n ^j . Injecting these inequalities in the right hand side of ( |56[ ) we 
conclude that 

2 



exp(e n ) > 



0n(2<5„) 



EL 

2 



4>n | 5n 
<5n <j>n 



> 1p 



1.36. 



where for the second inequality we used that x + x 1 > 2 for all x > 0. 
Let n, n + 1, ■ ■ ■ ,n + fc - 1 be consecutive integers in £+. Then 



fc ln(2p 2 



n+fe-l 

< E e * 



> i, 



In 



< In 



sin < 



sm </> n+fe 
sin(vr/2) 



+ 2 



cotan $ 



sin(2e) 



+ 2|lne| 
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where we used in the third line the hypotheses (i) and (iii). Regarding the last line, we have by 
concavity sin(2e) > (2/-7r)e > e, thus — lnsin(2e) < | lne|. The maximal number k of consecutive 
elements in E + is therefore bounded by c\\ lne|, where c\ := 3/ln(2p 2 ). 

The set E + can be arranged into at most #(E \ E + ) + 1 < cq + 1 series of consecutive 
elements, and by the previous argument they have length at most c\ \ lne|. Finally 



N = 2 + #(E \ E+) + #(E+) < 2 
which concludes the proof of this Lemma. 



co 



(c + l)ci| In el, 



□ 



Proposition 3.12. There exists a constant C such that the following holds. Let F E J, and let 

G be the norm defined by the diagonal matrix of entries k), k := Ak(F) + 1. 

Then there exists r < C(l + In k(F)) and 9\, ■ ■ ■ ,8 r £ M such that denoting tpi := <pg(' — 0%) 
one has 

min{(/9i, • • • , tp r } < tpp < max{(/9i, • • • , tp r } (57) 

Proof. We define $ e C ([0,2tt], [0,vr]) by *(t) := vr/2 - <p F (t). The difference $ - Id is non- 
increasing since tpp + Id is non-decreasing, see Proposition 3.6 For all t 6 [0, 27r] we have 
|¥>f(*)| < arccos(l/fc(-F)), using again Proposition 3.6, hence <5(i) > arcsin(l/K(.F)) > 1/k{F). 
Last for < t < t! < 2tt, one obtains using (40) 

ft' ,-t' 



\ cotan $ = / tan tpp < ]hk(F). 
Jt Jt 



Therefore $ satisfies the assumptions of Lemma 3.11, with e := min{-7r/6, 1/(2k(F))}. Let 
toi" ' >tN be the finite sequence given by this Lemma. We have N < C\ lne| < C'(l + \nn(F)) 
for some absolute constant C . Using (55) and the fact that $ — Id is non-decreasing, we obtain 
<fr(i) > e + (t — t n ) for all < n < N, and all t n < t < t n+ \, and therefore 

tp F (t) <ir/2-e- (t-t n ). 



We next observe, using Lemma 3.9 and with := arctan(l/K), that 



V?g(@) = arctan 



k — K 



> 



7T 



K — K 



7T 
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where we used successively that arctan(x) > ir/2 — 1/x for all x > 0, and n — k 1 > (Ak(F) + 
1) - 1 > &k(F) > 2/e. Therefore (p G (G + 1) > ir/2 - e - t, for all t > 0, using Q. Denoting 
#i := — for all < i < N, and (fi := v?g(- ~~ $i)> we conclude that 

(p F < max{93 ,--- ,<Pn}, 

on [0, 27r], hence also on M by 27r-periodicity. We have obtained one side of the announced 
inequality (57). 
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For the other side we define := ir/2 + (pp(—t), obtain to,- - ■ using Lemma 

with again N < (7(1 + In k(F)) for some absolute constant C. Setting 0i := — B + t{ anc 

<Pi := <Pg(' ~ iov < i < N, we obtain likewise (fp > min{<^i, • • • , This concludes the 

proof with r = N + N. □ 

We conclude in the next corollary the proof of main result of this paper: the average estimate 
of #(T(F S )), 6 G [0, 2-71"], for an arbitrary asymmetric norm F. 

Corollary 3.13. There exists a constant C such that for any asymmetric norm F on M 2 , one 
has 

\ #(T(F e ))al9 <C '(1 + In k(F)) 3 . (58) 
Jo 

Proof. In view of Lemmas |2.8| and |2.11[ we can assume that g £ Let k, G, r and 9i, ■ ■ ■ ,9 r 
be as in Proposition 3.12 Applying successively Proposition 3.8 and Corollary 3.10 we obtain 



r#{T{F e ))d9 < [ * #(T(G 6i+6 )) 

J0 l<i<r J ° 



dO 



Ki<r 



rC(l + 21iik) 



Recalling that r < C'(l + In k(F)) and k = 4k(F) + 1, see Proposition 3.12, we obtain the 
announced result. □ 



4 Implementation and Numerical results 

We compare our algorithm with two alternative solvers of the Escape Time problem, or Anisotropic 
Eikonal Equation, which enjoy a reputation of efficiency an simplicity in applications [lj: the 
Adaptive Gauss Siedel Iteration (AGS^J of Bornemann and Rasch |2j , and Fast Marching using 
the 8 point stencil (FM-8). Stencils are illustrated on Figure [2} Two more recent methods were 
also implemented: the Monotone Acceptance Ordered Upwind Method (MAOUM) of Alton and 
Mitchell [12] is tested when its memory usage allows it, and Fast Marching using Lattice Basis 
Reduction (FM-LBR) of the author [5] in the special case of Riemannian metrics. 

We considered four test cases: two involving (asymmetric) Finsler metrics, and two involving 
Riemannian metrics. Three of these tests are directly motivated by applications, including mo- 
tion planning, seismic imaging, and image segmentation, while the fourth one has the advantage 
of having an analytic solution, avoiding the recourse to a reference solution. Depending on the 
test, the metric anisotropy k{F) ranges from 4 to 400. Each algorithm was executed on each 

6 With stopping criterion tolerance parameter 10 -8 , as in [2]. 
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Figure 10: Level lines of the first test case [8]. L°° error (center) and averaged L 1 error (right), 
with respect to a reference solution, plotted as a function of CPU time, in this test case. Data 
obtained by running the algorithms at 100 different resolutions n x n, with n ranging from 61 
to 1201. Best viewed in color. FM-ASR: blue. FM-8: brown. AGSI: green. MAOUM: orange. 



See Figures 11 12 and 13 for the other tests. 



test case, at 100 different resolutions n x n, where n ranged from 61 to 1201 (odd values of n 
are preferred for symmetry reasons). We compare the algorithm's efficiency by representing, on 
Figures 10, 11 12 and 13 their L°° or averaged L 1 numerical error, with respect to an exact or 
reference solutiorj^J as a function of CPU time. We also discuss accuracy resolution-wise in the 
text: which resolution n is required to meet a prescribed L°° error bound ? 

In practical applications [I], the choice of the FM-8 versus the AGSI is typically regarded 
as a compromise in favor, respectively, of CPU time or numerical accuracy. Indeed the FM-8 
is a single pass solver with a small stencil, which thus completes in short and predictable CPU 
time, almost independent of the problem solved. For instance the FM-8 completes our four 
benchmarks in CPU tim^] 0.77s, 0.79s, 0.86s and 0.79s respectively, on a 601 x 601 grid, while 
the AGSI takea^] 8.61s, 285s, 15.1s and 123s. On the other hand the results produced by the 
AGSI are known to converge towards the viscosity solution of the continuous problem, as one 
refines the discretization grid, whereas convergence can only be guaranteed in limited cases for 
the FM-8. Indeed the acuteness condition (iii) in Definition 1 1 . 1 1 can be violated by the 8 point 
stencil, except for a Finsler metric J- such that k(J-) < y/2 (using (29)), a Riemannian metric 
such that ^(J 7 ) < y/2 + 1 (using Proposition 1.2 in |5J), or axis aligned anisotropy [11]. Our 



objective is to bring together the best of both worlds: our algorithm is fastF] (1.39s, 3.09s, 1.32s 
and 1.11s with the above settings) and universally convergent. 

The MAOUM [12] and the FM-LBR [5] are more recent algorithms than the AGSI or the 
FM-8, and are closer to FM-ASR from a theoretical point of view: they are universally conver- 
gent, inspired by Dijkstra's algorithm, and use static stencils assembled in a pre-processing step. 
However the MAOUM uses large isotropic stencils, see Figure [2j instead of smaller anisotropic 
stencils for the FM-ASR, resulting in a larger complexity and memory footprint. The FM-LBR 
mainly distinguishes itself from the FM-ASR by its domain of application: it is restricted to 
Riemannian metrics, but extends to higher dimension. 



Our first two tests involve asymmetric norms defined as the sum of an anisotropic euclidean 

7 Reference solutions computed on a 5001 x 5001 grid, using the AGSI in the first and third test cases, the 
FM-LBR in the last, and extended to the continuous domain via bilinear interpolation. 

8 A11 timings in seconds, obtained on a 2.4Ghz Core 2 Duo, using a single core. System memory: 8 GB. 
9 In the trivial case of a constant metric, equal to the euclidean norm, the AGSI takes 0.67s, on the same grid. 
10 CPU time for the FM-ASR, and the FM-LBR, includes stencil construction, which often accounts for 50%. 
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Figure 11: Level lines are circles in the the second test case. The shortest path joining reg to 
the origin, where eg := (cos 8, sin 6), is the spiral ip t-t (r — (p)ee- v , p £ [0,r] (red curve). 



norm and of a linear form. 

Proposition 4.1. Let M £ and let uj £ -ZR 2 6e such that (uj, Mu) < 1. T/ie map F : M 2 ^ 1R 
defined by 

F(u) := y/(u,Mu) - (uj, Mu) 

is an asymmetric norm, which unit ball {z; F(z) < 1} is an ellipse, not centered at the origin if 
uj ^ 0. Furthermore the dual asymmetric norm F* has the same form, with parameters M*,ui* 
defined by 

S:=l-lu,Mu), M* := ^ + J M , uj* := -M^uj/S. 



The minimization problem (12), appearing in the Hopf-Lax update operator, and the evaluation 
of the predicate "u,v form a F -acute angle", cost numerically 0(1) operations +, — , x, /, yf- 
among reals. 



Proof. Up to a linear change of variables (by Ma), we may assume that M = Id, and thus 
F{u) = \\u\\ — (uj,u), with \\u\\ < 1. The 1-Homogeneity and the Convexity of F are obvious. 
Furthermore F(u) > ||u||(l — > 0, with equality if and only if u = 0, which shows that 
F is proper, hence is an asymmetric norm. The boundary of the compact and convex set 
{z; F(z) < 1} is characterized by the quadratic equation ||ii|| 2 = (1 + (co,u)) 2 , which first degree 
term u i-> 2(uj, u) is non-zero if u) ^ 0. Hence this set is an ellipse, non-centered if ui ^ 0. 

Consider an arbitrary u £ H 2 \ {0}. Observing that 1/F*(u) = vnvn{F{v)] (u,v) = 1}, 
we find using the Khun- Tucker conditions for this constrained optimization problem, that the 
minimizer v satisfies v / \\v\\ — uj = Xu for some A £ H. Taking the scalar product of this equation 
with v we obtain A = 1/F*(u). On the other hand observing that ||cj + Xu\\ = 1, we obtain a 
quadratic equation which positive root is A. The announced expression of F*(u) follows. 

Since the norm F is differentiable, evaluating the predicate u u, v form an acute angle" is 
straightforward. It was observed from the very first works on fast marching methods |10j that in 



the special case ui = 0, the minimization problem (12) amounts to solving a quadratic equation. 



Choosing a non-zero oj instead is equivalent to subtracting (uj, u) to d u (resp. (uj, v) to d v ), thus 



the problem (12) has an equally simple solution. □ 



Our first test is a motion planning control problem, also discussed in [8j [2] . The Finsler metric 
is given by its dual: Fl(u) = [|u|| + (uj(z), u), where ui(x, y) = — 7sin(47rx) sin(47ry) and 7 = 0.9. 
The speed profile {u; F z (u)}, in the control theoretic interpretation, is the euclidean unit ball 
translated by uj(z), see Figure [l] (right): this could model a boat, able to move at unit speed 
on a still water, but caught in ocean current of speed uj(z). We compute the distance, i.e. the 



31 



minimal travel time, to the center of the square domain [—0.5, 0.5] 2 . See e.g. [5] for a discussion 
on shortest path extraction. The maximum anisotropy ratio, ^(F) = (1 + 7)/(l — 7) = 19, 
is not small, but anisotropy is pronounced only on a small region, where | sin(47rx) sin(47ry)| is 
close to 1. As a result, the FM-8 delivers excellent results in terms of averaged L 1 error, and 
best results in terms of L°° error for CPU times < Is, after what (presumed) non-convergence 
begins to show and the FM-ASR outperforms it. The FM-ASR is the best method among those 
which benefit from a convergence guarantee: for the prescribed tolerance 5 x 10~ 3 on the L°° 
error, the AGSI takes 11.3s, at resolution n = 661, while the FM-ASR takes 0.51s, at resolution 
n = 375, thus reducing CPU time by a factor 22. 

Our next test case is based on the following proposition. 

Proposition 4.2. Let g G C°(iR+,] — 1, 1[) be such that g(0) = 0. Let T be the Finsler metric 
defined by Tq(u) = \\u\\, u G M 2 , and for all z G M 2 \ {0} 

T z (u) = \\u\\- 9 -^-(z ± ,u), (59) 

where z 1 - denotes the rotation of z by 7r/2. Then the length D(z) of the shortest path joining z 
to the origin, solution of the eikonal equation ^ on ft := M 2 \ {0}, has the expression 

D{z)= I " ^l-g(r) 2 dr. (60) 



Proof. We consider z = roj, where r > and ||u;|| = 1, and we denote V(z) := g(r)w^ — 
y/1 -g(r) 2 uj. We have F z {u) + y/1 - g{r) 2 (u, u) = \\u\\ - (V(z),u) > 0, since \\V(z)\\ = 1, with 
equality if u is positively proportional to V(z). 

Let 7 G C ,1 ([0, 1],R 2 ) be such that 7(0) = z, 7(1) = 0. We may assume that 7(f) / for 
t < 1, up to eliminating a loop starting and ending at the origin at the end of the path 7. For 
t G [0, 1[ let r(t) := \\-y(t)\\ > and let u(t) := j(t)/r(t). Note that (co(t),j'(t)) = r'(t), for 
t G [0, 1[. We obtain using our lower estimate on F z {u) 



length(7):=/ F^ t) {i '(*)) dt > - y/\ - g{r{t)) 2 (oj(t), 7 '(t)> dt = / s/l - g(r{t)f r'(t) dt 



C ^/l-g(r{t)Y(u{t),i{t))dt= f 
Ji 

with equality if 7'(i) is positively proportional to V{^{t)) for all t £ [0, 1[ (a path of minimal 
length can therefore be obtained by solving an ordinary differential equation). Observing that 



the right hand side of (60) and of the last equation coincide, we obtain the announced result. □ 



For our second test case, we chose the Finsler metric T given by g(r) := r/Vl + r 2 in (59), 
in such way that D{z) = arcsinh(||,z||). The problem was discretized on a the square domain 
Q := [— To,ro] 2 , where ro = 10, which contains the disk B := {z; \\z\\ < tq}. Due to the 



allure of the paths of minimal length for the continuous problem, spirals, see Figure 11 the 
convergence of the discrete solution towards the continuous one can only be guaranteed within 
the disk B. Grid points that do not belong to B are thus rejected when computing errors. 
The maximum anisotropy ratio on the disk B is k(J-\b) = (ro + y/l + r^) 2 ~ 402, which is 
quite pronounced, and unsurprisingly the FM-8 non-convergence shows early. Due to the strong 
anisotropy, the MAOUM produced huge stencils, leading to a memory footprint incompatible 
with our equipment. Unlike other test cases, the AGSI produced here the most accurate results 
resolution-wise: for the prescribed tolerance 5 x 10~ 2 on the L°° error, the AGSI takes 119s, at 
resolution n = 435, while the FM-ASR takes 10.7s, at resolution n = 1069. Despite the higher 
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Figure 12: Test case inspired by seismic imaging (taken from [5], Figure 6, top left), with 
moderate anisotropy ^(J 7 ) = 4. Riemannian metric T z {u) := y/ (u, A4(z)u) , where A4(z) has 
eigenvalues 0.8 _2 ,0.2 -2 , the former associated to the eigenvector (1, (vr/2) cos(4-7rx)). Domain 
[-0.5, 0.5] 2 . Target L°° error bound 2 x 10 2 is met by the AGSI in CPU time 4.3s, at resolution 
n = 375, and by the FM-ASR in CPU time 0.18s (24 times less), at resolution n = 239. 




Figure 13: Test case inspired by tubular image segmentation, with strong anisotropy k(J-) = 
100, Figure 4 in pp. The Riemannian metric is equal to the euclidean norm, except on a band of 
width 1/100 along a spiraling curve V where it has eigenvalues 1/100 2 , 1, the former associated 
to the tangent vector to V. See [5] for details. Target L°° error bound 0.5 is met by the AGSI 
in CPU time 1015s, at resolution n = 1135, and by the FM-ASR in CPU time 0.054s (22700 
times less), at resolution n = 157. 



resolution, the FM-ASR strongly reduces CPU time needed to achieve a prescribed L° 
bound, here by a factor 11. 



Two more test cases, involving Riemannian metrics, are illustrated on Figures 12 and 13 



They are inspired by seismic imaging and medical image segmentation respectively, and were 
originally proposed in [8] and [I], see also [5]. The efficiency of the FM-ASR and of the FM-LBR 
[5] are comparable, and their superiority over alternative methods is here unquestionable. In 
the last test these two methods are in a class of their own, often reducing CPU time by four(!) 
orders of magnitude over their alternatives, for a target L°° error bound. 



Conclusion 

We introduced in this paper a variant of the fast marching algorithm, which applies to arbitrary 
Finsler metrics, on two dimensional domains discretized on a grid, and which is particularly 
efficient in the context of large anisotropies. Its complexity depends only (poly-)logarithmically 
on the anisotropy ratio k(J-), in an average sense over grid orientations, whereas earlier methods 
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had a linear or polynomial dependence in this parameter. Numerical experiments show a reduc- 
tion by an order of magnitude, or more, of the CPU time required to meet a target error bound. 
Future work will be devoted to the analysis of the accuracy of this algorithm, its extension to 
higher dimensions and triangulated domains, and its application to image analysis. 



Appendix : Proof of Proposition 1.3 



The optimization problem of interest ( 12 ) is the minimization of a continuous convex function 
on a compact interval, hence there exists at least a minimizer. Let G be the asymmetric norm 
defined by G(x,y) := F(xu + yv), for all (x,y) G K 2 . Let also D := (d u ,d v ) and 1 := (1,1). 
The problem (12) is equivalent to 



min{(w, D) + G(u); uj G R+, (uj, 1) = 1}. 



(61) 



The assumption that and 1 are not minimizers of the original problem (12), implies that the 
minimum (61) is not attained when to is equal to e x '■= (1,0) or e y := (0,1). We denote by 
uj a minimizer of (61), and observe that both components of to are positive. The Kuhn- Tucker 
relations, for this constrained optimization problem, state that there exists a scalar A G H, the 
Lagrange multiplier, and an element V G dG(oj) such that 



D + V = \1. 



(62) 



We denoted by dG(to) the sub-gradient of the convex function G at the point uj; if G is differ- 
entiable at uj, then dG(uj) = {VG(w)}. Since G is 1-homogeneous, we have (to, V) = G(uj), by 
Euler's homogeneous function theorem. Taking the scalar product of (62) with uj we obtain 

A = \(uj, 1) = (uj, D) + (to, V) = (uj, D) + G(uj) = d w . 



Injecting this relation in (62), we obtain V = d w l — D. In order to conclude the proof, we need 
to show that both d w — d u = (e x , V) and d w — d v = (e y , V) are positive. Since the minimum 
(61 ) is not attained for uj = e x , we have 

d u + G(e x ) >d w = X(e x , 1) = (e x , D) + (e x , V) = d u + (e x , V), (63) 

thus G(e x ) > (e x , V). On the other hand, denoting by (a, f3) the components of uj, and recalling 
that they are positive, we obtain 

a(e x , V) + f3(e y , V) = (uj, V) = G(uj) = G(ae x + f3e y ) > aG(e x ), 

where for the last inequality we used that e x , e y form a G-acute angle, since u, v form a F-acute 
angle. Finally 

G(e x ) > (e x , V) > G(e x ) - ((3/a)(e y , V), 
hence (e y , V) > 0. Likewise (e x , V) > 0, which concludes the proof. 
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